Non-local pair correlations in the ID Bose gas at finite temperature 
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The behavior of the spatial two-particle correlation function is surveyed in detail for a uniform ID 
Bose gas with repulsive contact interactions at finite temperatures. Both long-, medium-, and short- 
range effects are investigated. The results span the entire range of physical regimes, from ideal gas, 
to strongly interacting, and from zero temperature to high temperature. We present perturbative 
analytic methods, available at strong and weak coupling, and first-principle numerical results using 
imaginary time simulations with the gauge- P representation in regimes where perturbative methods 
are invalid. Nontrivial effects are observed from the interplay of thermally induced bunching behavior 
versus interaction induced antibunching. 
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INTRODUCTION 



The study of two-body correlations has a long history 
dating back to the 1956 experiment of Hanbury Brown 
and Twiss (HBT) fl. The HBT experiment set out to 
measure the intensity of light coming from a distant star, 
at two nearby points in space. The fluctuations in the in- 
tensities were shown to be strongly correlated in spite of 
the thermal nature of the source. In more recent times, 
experimental progress in the field of ultra-cold atomic 
gases has provided the opportunity to examine similar 
correlations in systems of cold atoms (as opposed to pho- 
tonic systems) . The large thermal de Broglie wavelength 
in a cold gas means the correlations occur on length scales 
large enough to be resolved using current detectors. A pi- 
oneering experiment of this kind involving a cloud of cold 
Neon atoms, was carried out by Yasuda and Shimizu [2] 
as early as 1996. A more comprehensive study was un- 
dertaken during 2005 - 2007 in Refs. @, 0], where the 
two particle bunching phenomena associated with Bose 
enhancement (when metastable 4 He* atoms were used) 
was juxtaposed with the antibunching behavior present 
in a system of fermions (when 3 He* atoms were used). 
In all of the above cases the measured correlations were 
completely described by the statistical exchange interac- 
tion between particles in an ideal gas. 

The behavior of strongly interacting systems poses 
some of the most difficult questions confronting current 
theoretical studies in many-body physics. In this paper 
we discuss how our simple understanding of two-body 



correlations in an ideal gas can be radically altered in 
the presence of interactions. To demonstrate this we cal- 
culate the normalized pair correlation function 



g V)(r) = (¥{Q)W(r)<S>{r)<S!(Q))/n 



(1) 
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in a homogeneous repulsive one-dimensional (ID) Bose 
gas 0, Q at finite temperature over a wide range of in- 
teraction strengths. In Eq. JT]), ^(x) is the field operator, 
and n = (x)if* (x)) is the linear ID density. Physically, 
gr( 2 )(r) quantifies the conditional probability of detecting 
a particle at position r, given that a particle has been 
detected at the origin. Theoretically the ID Bose gas 
model with <5-function interaction is one of the simplest 
paradigms we have of a strongly interacting quantum 
fluid, owing to its exact integrability [E IE 0, H IE EE • I n 
the limit of an infinitely strong interaction it corresponds 
to a gas of impenetrable (hard-core) Bosons treated first 
in Ref. [ll|. It also holds relevance as an experimen- 
tallyaccessible system [H, El EI EI EE E3, EE El, HE 
mM HI HI EE HE- Opposite from 2D and 3D, the 
strongly interacting limit of a ID system is achieved in 
the low density regime. In this regime the wave function 
of the particles is strongly correlated and prevents them 
from being close to each other, which results in dramatic 
suppression of 3-body losses. This allows for the stable 
creation of strongly interacting ID Bose gases. 

There has been a substantial amount of previous the- 
ory on correlations of the ID Bose gas model. The Lut- 
tinger liquid approach provides a method of calculating 
the long-range asymptotic behavior in the decay of non- 
local correlations |9|, LlC|| - Local second- and third-order 
correlations in the ho mog eneous system have been calcu- 
lated in Refs. [IE HE HE HE HH j extensions to inhomo- 
geneous systems using the local density approximation 
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(LDA) are given in Ref. [32|. Numerical calculations at 
specific values of interaction strength have been carried 
out at T = [33| and at finite temperature j34[ . Similar 
nonlocal quantities have been calculated for the T = 
37,[H,[3§], and for finite temper- 
34 1 and in the strong interaction 
contain recent re- 



ground state [331.135 
ature both numerically 
limit [13]. Refs. @, ©, Efl E3, S3, SI 

views of the physics of the ID Bose gas problem. 

The focus of the present paper is the nonlocal corre- 
lation function at arbitrary interparticle separations r; 
we give the details of analytic derivations of the results 
discussed in a recent Letter [3| and complement them 
with exact numerical calculations using the stochastic 
gauge-P method of Ref. H SI, S3, [H]. Experi- 
mental proposals to measure nonlocal spatial correlations 
between the atoms in a ID Bose gas have been discussed 
in Ref. 0521]. 

The structure of this paper is as follows. In section ITT1 
we give a brief review of the physics of a ID Bose gas, em- 
phasizing the important parameters which determine the 
phase diagram. In section [Hi] we outline the details in- 
volved in the application of the (imaginary time) gauge-P 
phase space method to the ID Bose gas. The more tech- 
nical details are placed in appendix [Aj This method is 
capable of obtaining numerical results in the cross-over 
regions of the phase diagram, where analytic results are 
not available. In sections IIV} FVl and IVII we present the 
results of calculating g^ 2 '(r) in the nearly ideal gas limit, 
the weakly interacting limit, and the strongly interacting 
limit respectively. The results are obtained from numer- 
ical calculations and analytic perturbation expansions. 
We describe the details of our perturbation expansion 
in each respective section. In section IVIII we analyze, 
in detail, the nature of the crossover into the fermion- 
ized Tonks gas regime. Section [Villi discusses the limita- 
tions of the numerical method. In section HXl we give an 
overview and draw conclusions. 



II. THE INTERACTING BOSE GAS IN ID 

We are considering a homogeneous system of N identi- 
cal bosons in a ID box of length L with periodic bound- 
ary conditions @, @|. We include two-body interactions 
in the form of a repulsive delta-function potential. The 
second-quantized Hamiltonian of the system is given by 



H = — [ dxd x ¥d x ^> + - [ dx¥¥M, 
2m J 2 J 



(2) 



where m is the mass and g > is the coupling con- 
stant that can be expressed via the 3D s-wave scattering 
length a as g ~ 2h 2 a/{mt 2 1 ) = 2huj±a [50]. Here, we 
have assumed that the atoms are transversely confined 
by a tight harmonic trap with frequency u>± and that 
a is much smaller t han the transverse harmonic oscil- 
lator length Zj_ = y^h/mui^. The ID regime is realized 
when the transverse excitation energy hui± is much larger 
than both the thermal energy T (with fc^ = 1) and the 



chemical potential /i [32j, |51| . A uniform system in the 
thermodynamic limit (N, L — > 00, while the ID density 
n = N/L remains constant) is completely characterized 
@; 0] by two parameters: the dimensionless interaction 
strength 



7 



mg 
h 2 n 



and the reduced temperature 

t = T/T d , 



(3) 



(4) 



where Td = h 2 n 2 /(2m) is the temperature of quantum 
degeneracy in units of energy (3(J • 

The interplay between these two parameters dic- 
tates the dominating behavior in six physically different 
regimes. Briefly, these regimes are: 

• Nearly ideal gas regime, where the temperature 
always dominates over the interaction strength. 
This regime splits into two subregimes defined by 
t « 1 or t > 1. In both cases one must have 
7 -C min {r 2 , y^}- 

• Weakly interacting regime, where both the interac- 
tion strength and the temperature are small, but 
r 2 < 7 « 1. This regime realizes the well known 
quasi-condensate phase. Fluctuations occur due to 
either vacuum or thermal fluctuations, which de- 
fines two further subregimes, with r < 7 or r > 7, 
respectively. 

• Strongly interacting regime, where the interaction 
strength is large and dominates over temperature 
induced effects. This can occur at high and low 
temperatures, again defining two subregimes with 
t < 1 or r > 1. 

The basic understanding of the competition between in- 
teraction induced effects and thermally induced effects 
was outlined in Ref. 

Although the model is integrable via the Bethe ansatz, 
the cumbersome nature of the eigenstates [H3| inhibits 
the direct calculation of the nonlocal two-body correla- 
tion function. We therefore use numerical integration 
in a phase-space representation, together with perturba- 
tion theory in each of the six regimes. The standard 
Bogoliubov procedure, applied to Eq. (J2j) is appropriate 
in the case of the weakly interacting regime (see section 
IV|) . Perturbation theory in the strongly interacting and 
nearly ideal gas regimes is done using the path integral 
formalism (see sections llV Al and IVII respectively) . 



III. NUMERICAL STOCHASTIC GAUGE 
CALCULATIONS 

A. Gauge-P distribution 

To evaluate correlations away from the regimes of ap- 
plicability of the analytic approximations, we use the 
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gauge-P phase-space method to generate a stochastic 
evolution from the simple T — > oo limit (where inter- 
actions are negligible) down to lower temperatures. This 
method gives results that correspond exactly to the full 
quantum mechanics using the Hamiltonian ([2]) as the 
number of averaged realizations (S) goes to infinity. The 
gauge- P method has been described in [H, EE 53], and 
is covered in greatest detail in [i^], while an initial ap- 
plication to the ID Bose gas was presented in [34)]. Be- 
low we give a summary of the derivation for this system, 
and present the basic calculation procedure. Some of the 
more technical details are given in Appendix [A] 

We consider a grand canonical ensemble with mean 
density n, Hamiltonian ([2]) and inverse temperature given 
by (5 = 1/ksT. When the Hamiltonian commutes with 
the number operator N = J dx^ (x)^> (x) , as is the case 
here, the unnormalized density matrix at temperature T 
is given by 



■2 - p W)N-H]p 
Pu — e , 



(5) 



where p{P) is the chemical potential. In this formulation, 
p, can in principle be chosen at will as any desired function 
of temperature, thus indirectly determining the density 
n(T). In the Schrodinger picture the density matrix is 
equivalently defined by an "imaginary time" master-like 
equation 



dMP) 
dp 



p e {p)N-H\ p u (P) 
p e (P)N-H , p u (P) 



(C) 



([6]), and then in a second step, sampling this distribu- 
tion stochastically and evolving the samples with a diffu- 
sive random walk that is equivalent to the Fokker-Planck 
equation. The general approach is described in [53L l54l|. 
The price that is paid for tractable calculations is a loss 
of precision that comes about due to the finite sample 
size S. Fortunately this uncertainty can be readily es- 
timated using the Central Limit theorem and scales as 

Vs. 

We utilize the normalized off-diagonal coherent state 
expansion of the positive-P distribution [H3| because the 
number of variables required to describe a sample is linear 
in the number of spatial points (tractability) and because 
it describes all quantum states with a non-negative real 
distribution. However, for this investigation two addi- 
tional elements are needed. Firstly, the evolution ([6]) does 
not preserve the trace, so an additional weight variable in 
the expansion is needed to keep track of this. Secondly, 
the evolution equations for the samples given by a bare 
weighted positive-P treatment are unstable and can lead 
to systematically bad sampling [H3|. The complex part of 
the weight variable allows us to remove these instabilities 
using a stochastic gauge as described in (3~5 . |45| . 

In practice, the first step is to discretize space into 
M equally spaced points in a box of length L with peri- 
odic boundary conditions, on which the fields are defined. 
There is a lattice spacing of Ax = L/M per point. One 
must make sure that the lattice is fine enough and long 
enough to encompass all relevant detail. In practice we 
check this by increasing L and, separately, M until no 
further change in the results is seen. Having this equiv- 
alent lattice, one can expand the density matrix p u as 



and a simple initial (i.e. T — > oo) condition 



p u = I G(v)A(v) d 4M+2 v, 



(9) 



Pu(0) 



-XN 



(7) 



with A = — lim^o [Pp-(P)} and /? playing a similar role 
to time in the Schrodinger equation for time evolution, 
apart from a factor of i (hence the name). The sec- 
ond line of ([6]) follows from the restricted set of density 
matrices described by the grand canonical ensemble ((!]), 
where logp u commutes with p u . Note that p e (P) is a 
temperature-dependent "effective" chemical potential 



Me 



dp ' 



(8) 



that is not necessarily equal to p. The initial condition 
§2§ can then be evolved according to Eq. (J6J to obtain 
the equilibrium state at lower temperatures P > 0. How- 
ever, in the density matrix form, this naturally becomes 
intractable for more than a few particles. 

Phase-space methods such as the gauge-P distribution 
used here reduce the computational resources needed to a 
manageable number. This is done by deriving a Fokker- 
Planck equation for a distribution of phase-space vari- 
ables that is equivalent to the full quantum mechanics 



with a positive [45] distribution G(v) of the set of 2M + 1 
complex phase-space variables, 



, «M, Oil, 



r>«}> 



(10) 



v = {ai, . 
that describe an operator basis 

A(0) = n®^ 1 IK)((a+ni e-EJW * (11) 

composed of unnormalized (Bargmann) coherent states 
= exp OjV Ax W(xj) |0) at the j-th point at lo- 
cation xj = (j — l)Ax and a global weight Q. 

The initial condition ((7|) corresponds to the distribu- 
tion 



M 



+ exp(-|a i | 2 /n x ) 



G (v)^S^n-l)l[5 2 (oy-Cat) 

7=1 

(12) 

where n x = l/(e A - 1) = N/M is the mean number of 
atoms (N = (N)) per spatial point in the initial P = 
state. We see that, at least initially, a + = (a)* are 
complex conjugates. 
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B. Fokker-Planck Equation 

To generate the Fokker-Planck equation (FPE) for 
G(v) corresponding to the master equation ^ we use 
the following differential identities for the basis opera- 
tors 



/As*(ij)A = Oj-A 
/ Ax$ t (x- J )A : 



3 " I1 J 

a+- 



\fKxMf(x 3 ) = a+A, 



d 



dc 



A. 



A. 



(13a) 
(13b) 
(13c) 
(13d) 



These convert quantities involving the operators ^, ^ 
and p u to ones involving only A and their derivatives. 

In what follows it will be convenient to label the a and 
a + variables as 



which is initially the number of particles at the j-th site, 
and an effective complex-variable Gibbs factor K corre- 



sponding to Tr 



(H-fi e N)A /Tr 



A 



K(v) = 



h 2 (Va+) (V aj 
2m 



9_Nf 
2Ax 



(16) 



Here Voj is the discretized analogue of the gradient of a 
complex field a{x) that satisfies ct(xj) = ay. 

To obtain a FPE equation for G(v) we proceed as fol- 
lows. Firstly, we can make use of the additional "gauge" 
identity that follows trivially from Eq. ifTTj) . 



(17) 



» 



j ■ 



if v = 1, 
if v = 2. 



Using i|13p on ([6]) one obtains 
dG(if)_£ dAM+2 _ t _ 



d/3 



v = - G(v) 



4Ax 



E(«r) 2 - 



with 



K{v) 





\ — + 




) 9 




j daj 


V da 3 , 


iDaf 



N, 



(14) 



Ad 4M+2 v, 



(15) 



to convert K{v)A = K(v)Q-J^A on the first line of Eq. 
(fl"4"l) . This step is necessary in order to obtain an equa- 
tion of a form that can later be sampled with a diffusive 
process. Secondly, we integrate by parts to obtain differ- 
entials of G rather than A. Thirdly, if the distribution 



G is well bounded as la,- 



mi 



oo, we can discard 



the boundary terms. As it turns out (see appendix I A 1[) . 
this is not fully justified for the equation (fl"4)l . and the 
boundary behavior will need to be improved with the 
help of a stochastic gauge as described originally in [45| . 
However, for demonstrative purposes let us proceed on 
for now, and return to remedy the problem below in 
Sec. IIIIDI Lastly, having now an equation of the form 
J Ax [Differential operator] G(v) dv = 0, one solution is 
certainly [Differential operator]G(w) = 0, which is the 
following FPE: 



an 



VLK{v) 



d_ 

df3 



E 

3,v 



(«r) 2 



i d ( h 2 (y 2 a 



2 daP 



j > 



2m 



(«0 



Ax 



G(v). (18) 



C. Equivalent diffusion 



to the following set of stochastic differential equations 



A diffusive random walk that corresponds to the 
Fokker-Planck equation lfT8|) is found by replacing the 
analytic derivatives with appropriate derivatives of the 
real and imaginary parts of oq 1 ' [13, [f54[. This results in 
a diffusion matrix in the phase-space variables v with no 
negative eigenvalues. In the Ito calculus this is equivalent 



da 



dd 



dn 

df3 



2 I Me 



0) 



h 2 v 2 

2m 



9NA a M 
Ax J 3 



9 A») 



(19) 
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We do not use diffusion gauges [47| here and decompose 
the diffusion matrix in the most straightforward fashion. 
Here, the Q (fl) are real, delta-correlated, independent 
white Gaussian noise fields that satisfy the stochastic av- 
erages 



O as follows: 



= o, 

= 8ij8 vv >5(fi - 



(20a) 
(20b) 



In practice, at each time step separated from the subse- 
quent by an interval A/3, one generates M independent 
real Gaussian random variables of variance 1 / A/3 for each 

A») 
S ' 

Equations (fl9|) can be intuitively interpreted by noting 

that the equation for the amplitudes at each point is 
a Gross-Pitaevskii equation in imaginary time, with some 
extra noises that emulate the wandering of trajectories in 
a path integral formulation around the mean field solu- 
tion given by the deterministic part. A different wander 
for different v. The weight evolution of f2 generates the 
Gibbs factors of the grand canonical ensemble. 



D. Final equations 

A straightforward application of the diffusion equa- 
tions lfT9|) is foiled by the presence of an instability in 



the da.j /df3 equations. We use a stochastic gauge to 
remove this instability, in a manner described in |47tl4gj|. 
with the details given in Appendix IA 11 The final Ito 
stochastic equations of the samples are 



da 



d/3 




ilmNj) 



(21) 



dp 



-K{v) - i 



2Ax 



Some technical details regarding integration procedure, 
importance sampling, and choice of /Lt e (/3) are given in 
Appendix El Attention to these issues can speed up the 
calculations and reduce sampling errors by orders of mag- 
nitude. 



E 



O 



Tr 


Opu 


j G(v)Tr 


OA(v) 


dv 


Tr [p u ] 


JG(v)Tt 


~A(v) 


dv 





OA(v) 




d,v 




~A(vj 


\ Re(ft) 5 

/ 5 



(22) 



where (• • • )s denotes a stochastic average over the sam- 
ples, and T is an appropriate function of the phase-space 
variables v. The last line follows from properties of the 
operator basis A, and because the trace of p u and of ex- 
pectation values are real. 

The identities (fl3|) can be used to readily evaluate T 

since Tr A = Q. In particular, 



Re((A^)), 
Ax Re (fi) < 



(23) 



m^m^))= ^^f , (24) 



which explains the relationship between Nj and the par- 
ticle number at the j-th site. For the uniform system 
considered here, it is efficient to average the quantities 
over the entire lattice, so that e.g. 

L ( J &(x)&(x + r)$(x + r)$(x) dx^ 

9 {2 Hr) = ^ -a L - (25) 

U *t(a;)*(a;) dx\ 

Uncertainty is estimated as follows: We separate the S 
realizations into B bins, such that B 3> 1 and S/B ^> 1. 
One calculates an estimate for the expectation value of 
an observable in each bin independently (let us denote 
Oi as the estimate obtained from the zth bin). The best 
estimate for the expectation value of the observable is 
obviously (Oi)g. The one-sigma uncertainty in this es- 
timate is obtained from the Central Limit theorem and 



AO 



-n_J{0~) B -{0)% 



B 



(26) 



IV. 



NEARLY IDEAL GAS REGIME 

[7 < min{r 2 , ^t}] 



E. Evaluating observables 

Given S realizations of the variable sets v, using fresh 
initial samples and noises Q (/3) each time, one gener- 
ates an estimate of the expectation value of an observable 



We now present the perturbation theory results for the 
decoherent regime of a ID Bose gas [13], where both the 
density and phase fluctuations are large and the local pair 
correlation g^ 2 \0) is always close to the result for non- 
interacting bosons, </ 2 )(0) = 2. Depending on the value 
of the temperature parameter t, we further distinguish 
two sub-regimes: decoherent classical (DC) regime for 



t ^> 1 and decoherent quantum (DQ) regime for temper- 
atures well below quantum degeneracy, r < 1. Both can 
be treated using perturbation theory with respect to the 
coupling constant g around the ideal Bose gas, for which 
the nonlocal pair correlation function has been studied 
in Ref. Here, we extend these results to account for 
the first-order perturbative terms. 



A. Perturbation theory in 7 

The correlations of a ID Bose gas are governed by the 
action 

S[^*^] = J da Jdr [y*d a ^- W(¥*,#)], (27) 

written in terms of a space and imaginary time dependent 
c-number fields \&(x, a) in the Feynman path integral for- 
malism. Here a is the imaginary time and (3 = 1/fcsT 
is the maximum, corresponding to the inverse tempera- 
ture. The Hamiltonian density Ji is obtained from ([2]) by 
replacing the operators with the c-number fields. Using 
action lj27j) . the pair correlation function is given by 

g( 2 \r) = A- [ e- 5 [***l**(0)**(r)*(r)^(0). 
n Z J 

(28) 

where Z = JV^*^ e s ^ *1 is the partition function. 
In Eq. ([28]) and below, we use the notation that fields 
with imaginary time dependence omitted act at a = 0, 
i.e. ty(r) = ^(r, 0). Expanding the action pT)) in powers 
of g, we obtain up to the first order 

9 (2 Hr) =g£L(r) -^f da J dr ' 

x *(r',cr)*(r',cr)**(0)1'*(r)*(r)*(0)), (29) 

where c^liM = 1 + G(r, 0~)G(-r, CT)/n 2 is the ideal 
Bose gas result following from Wick's theorem. Note that 
since the expansion above is formally in powers of g, the 
final result can always be expressed in powers of 7 as 
7 oc g. The average in Eq. (|29|) is evaluated using Wick's 
theorem 156] 



A fl W(r) = g (2) {r) - 9$L(r) =~\ da I dr' (30) 

n Jo J 

x G(r',a)G(r - r' , -a)G(r' - r, a)G(-r',-a), 
with the Green's function 



2g 



G(r,a) = -(¥(0,0)** (7, a)) 



Akr—ihuj n (7 



f3L ^— ' ihuj n — h 2 k 2 /2m + /./ 



(31) 



The uj n ((3) are the Matsubara frequencies and the imagi- 
nary time a runs between and (3. The Green's function 
is periodic in the case of bosons and anti-periodic in the 



case of fermions. Thus it can be Fourier transformed with 
u) n = 2nn/P (bosons) or ui n = ir(2n + l)/(3 (fermions). 
The discrete sum over k becomes an integral in thermo- 
dynamic limit. 

In terms of a Green's function Gk(a) that is Fourier 
transformed with respect to the spatial coordinates, 
Ag( 2 >(r) can be brought to the form 

^ {2 \r)=- 2 ^j\j ^-e ikr T(k,a)r(k,-a), (32) 
where 



F(k,a) = — / dp G p+k (a)G p (-a), (33) 



and 



-n fc (/3)e-^ 2fe2 / 2 ™^), a<0, 
Gk ^>-^ -[l + n fe (/3)] e -^ 2 /2 m -M) 5 a>0j ^ 



with 



n fc (/3) 



(35) 



e (h 2 k 2 /2m-p.)l3 _ I 

being the standard bosonic occupation numbers. 

B. Decoherent classical regime 

For temperatures above quantum degeneracy, r ^> 1, 
the chemical potential is large and negative, so the 
bosonic occupation numbers are small, nfc(/3) <C 1, 
and can be approximated by the Boltzmann distribu- 
tion, rife(/5) ~ e -(h k /2m-iJ,)p Accordingly, the function 
Gk{a) in Eq. ([34]) becomes a Gaussian 



G k (a) 



cxp[-(ft 2 fc 2 /2rn- fi)(a + 0)], a < 0, 
-cxp{-(h 2 k 2 /2m~ fi)a} 7 a > 0, 



and Eq. 11331) is integrated to yield 



r(fc, a) = T(k, —a) = ne 



-a{0-u)h 2 k 2 /2m/3 



(36) 



(37) 



Here the mean density at a given temperature and chem- 
ical potential is determined from n = ^ J dk Gfc(0 _ ) = 
v / m/(2Trh 2 (3) e 13 ^. Using Eq. (021), the correction ((32j) to 
the pair correlation function is found as (see Appendix 



AgM(r) 




(38) 



where erfc(a;) is the complimentary error function. 

Together with g^Li( r ) = 1 + cxp[-™ 2 r 2 /2] (t > 1), 
this gives the following result for the pair correlation 
function in the DC regime (r 3> max{l,7 2 }): 



/ 2 )(r) = l + e-^/ A -) 2 



2tt 7 2 



erfc 



A-, 



(39) 
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FIG. 1: Nonlocal pair correlation g^ 2 \r) in the nearly ideal 
gas regime: (a) decoherent classical regime, r ^> max{l,7 2 }, 
Eq. ll39l) . with r in units of the thermal de Broglie wavelength 
At = \/4:TY/(rn 2 ); (b) decoherent quantum regime, y^y <C 
t <C 1, Eq. (1451) . with r in units of the phase coherence 
length lif, = 2/nr. 



This is written in terms of the thermal de Broglie wave- 
length 



A T = 



2nh 2 
2m,T 



in 



(40) 



a quantity that will appear repeate dly in what follows. 
At r = we have </ 2 '(0) = 2 — y^phnjr in agreement 
with Ref. [30| . In the non- interacting limit (7 = 0) we 
recover the well-known result for the classical ideal gas 
[H^ | characterized by Gaussian decay with a correlation 
length At- For 7 > we observe [see Fig. fflja)] the 
emergence of anomalous behavior, with a global maxi- 
mum g( 2 \r max ) = g( 2 \0) + 2j 2 /t at nonzero interpar- 
ticle separation nr max = 27/r <C 1. This corresponds 
to the emergence of antibunching, </ 2 ^(0) < g^ (r max ), 
due to repulsive interactions. As 7 is increased further, 
there is a continuous transition from the DC regime to 
the regime of high-temperature "fermionization" (see Sec. 
IVIB|) . with g^{0) reducing further and the maximum 
moving to larger distances. 



C. Decoherent quantum regime 

For temperatures below quantum degeneracy, with 
€ t « 1, only ui n = contributes to the Green's 
function 

G k (a)^-T[h 2 k 2 /(2m) + \ f i\}-\ (41) 

which gives the relation b etween the density and the 
chemical potential n = Tyjm/(2h 2 \[i\), [i = — Per- 
forming the Fourier transform of Eq. f4Tj) one obtains the 
one-particle density matrix for the ideal gas 

= <* f (0)*(r))/n = exp(-r/W, (42) 

which characterizes the decay of phase coherence over a 
length scale given by 

h 2 2 

h = n i r = — . (43) 
2m|/i| nr 

and also determines the second-order correlation function 
for the ideal gas 

|2 _ 1 , „-2r/l$ 



ffSLw = i+iflSLwi a = i+e- 



(44) 



The one-particle Greens function, Eq. (|4T|) . together 
with Eq. d3nj leads to T(k, a) = 4n%/(fc 2 Z 2 +4). Insert- 
ing it into Eq. lf32j) we obtain (see Appendix! 



correc- 



(2, 

tions to 5i,jeai( r )' leading to the following result for the 
pair correlation function in the DQ regime 

47 A 2r 
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(2 



) (r) = l + 



1 - 



1 



l<j> 

(0 



(45) 



2 - 4 7 /r 2 



This has the maximum value g K 

agreement with the result of Ref. [30j. For 7 = the 
correlations decay exponentially with the characteristic 
correlation length of half a phase coherence length de- 
scribing the long-wavelength phase fluctuations. 

An interesting feature in this regime is the apparent 
prediction of weak antibunching at a distance as seen in 
Fig.[T](b), with g' 2 '(r m ; n ) < 1. The strongest antibunch- 
ing in expression (|45|) occurs at rar m i n = r/47 3> 1, or 
J'min = l(f>T 2 /A-y ^> 1$, and dips below unity by an amount 
(47/T 2 ) cxp(— t 2 /47) <C 1. However, there is ambiguity 
regarding its existence: One should note that the dip 
below unity is very small in the region of uncontested 
validity of Eq. (|45f where r/y/y ^> 1, and only becomes 
appreciable around r < which is in the crossover 

region into the quasi-condensate (see Sec. [V]). Whether 
such anomalous antibunching survives higher order cor- 
rections in the small parameter y/j/r remains to be seen. 
Our numerical calculations to date have not been able to 
access a regime of small enough t/j/t to confirm or deny 
its existence. 

The numerical examples shown in Fig. [2] are for 
■y/7/V ~ 0.24 and y/y/r ~ 0.77, and show a thermal 
bunching peak with a typical Gaussian shape at the 
shortest range of A T , with A T <C 1$. At longer ranges, 
phase coherence dominates this and leads to exponential 
decay on the length scale 1$, in agreement with Eq. (|45|) . 
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FIG. 2: Approach of the pair correlation function to the 
ideal gas solution (shown dashed) in the decoherent quan- 
tum regime at r = 0.1, wit h r in u nits of the thermal de 
Broglie wavelength, At = \J 4-k/tii 2 . The thickness of the 
solid lines (numerical results) comes from the superimposed 
lcr error bars which are below resolution. 



D. Quantum/classical transition 

The transition from the quantum to the classical deco- 
herent gas was investigated using the gauge- P numerical 
method. The behavior is shown in Figs.EHU 

With rising temperature, still below degeneracy, one 
first finds a rounding-off of the exponential behavior at 
short ranges of a fraction of At, as seen in Fig. [21 There is 
also a global lowering of (r) with 7. It should be noted 
that the parameters for the numerical results shown in 
Fig. [2] are not deep in the regime where 1(45]) applies ac- 
curately, and the lowering of the tails with 7 is weaker 
here, than predicted by that limiting expression. 

Considering variation with T, as temperature ap- 
proaches, and then exceeds Gaussian thermal- like 
behavior appears first at short ranges, progressively tak- 
ing over an ever larger part of g^ 2 '(r) as temperature is 
raised. This is seen in Fig. [3j The exponential tails can 
persist at ranges r > At/V2tt well into the high tem- 
perature regime when 7 is small, as seen in Fig. [3jb) for 
t = 3 and even r = 10. 

There are three scenarios that can typically be con- 
trolled in ultracold gas experiments: (i) varying the ab- 
solute temperature changes r but not 7, as in Fig.[3j (ii) 
varying the coupling strength via a Feshbach resonance 
or varying the width of the trapping potential affects 7 
but not r, as considered in Section fVIII and Fig. O and 
(iii) varying the linear density gives changes in both 7 
and t, while keeping the quantity y/\/r constant. No- 
tably, this is the parameter that appears in the analytic 
expressions for both decoherent regimes, Eqs. (|45| and 

Figure[Hshows the behavior under scenario (iii), where 
increasing r corresponds to decreasing density of the 




FIG. 3: Exact behavior of g^ 2 \r), with r in units of At, in the 
nearly ideal gas regime with 7 = 0.001 and varying r around 
the quantum/classical crossover. In panel (b), the derivative 
/ = d[ln(g (2 \r) — l)]/dr shows a clear distinction between 
exponential decay (when / is constant) and Gaussian thermal- 
like behavior when / is linear. The triple lines indicate the 
numerical curves together with la error bars which are mostly 
below resolution. 



y/x 1/2 = 0.03 
1=1,10,100,1000 




gas. As expected, g^(0) tends to a constant value 
7a/27t/t 7^ 2 with r — > 00 predicted by Eq. 



ff (2) (0) = 

(|39)l . Interestingly, the crossover is quite broad under 



FIG. 4: Approach to the classical decoherent gas solution 
(shown dashed), Eq. (J39j) , for finite but small interaction 
with y/y/r = 0.03, which corresponds to a variation of density 
while keeping the coupling g and T constant. Here </ 2 ^(0) — > 
1.925 in the r — > 00 or equivalently n — > limit. Triple 
solid lines are the numerical results, with la error bars below 
resolution. 
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FIG. 5: Behavior of g {2) (r) in the crossover region between 
decoherent classical and quantum gas at r = 1. Values of 7 
shown are 0.001, 0.003, 0.01, 0.03, 0.06, 0.1 and 0.2 as the 
curves for g^ 2 \r) descend. 



changing density, with departures from the decoherent 
classical result still visible at r ~ 100. 

Finally, in the middle of the crossover region at r = 1 , 
7 « 1, there is the smooth and quite broad transition 
from low values of 7 to 7 ~ 0(1) that is shown in Fig. [H 
The situation of a short-range Gaussian with standard 
deviation ~ and exponential tails with length 

scale l$/2 that was seen in Fig.[3]morphs into an anoma- 
lous form with a local maximum that is similar to the 
high temperature fermionization behavior described be- 
low in Sections IVII and IVII1 



V. WEAKLY INTERACTING 
QUASI-CONDENSATE REGIME [r 2 < 7 < 1] 

In the regime of weak interactions and low tempera- 
ture (or Gross-Pitaevskii regime) with 7 <C 1 we rely on 
the fact that the equilibrium state of the gas is that of a 
quasi-condensate [Ha . [59| . In this regime the density fluc- 
tuations are suppressed while the phase still fluctuates. 
The pair correlation function is close to one and the de- 
viations can be calculated using the Bogoliubov theory. 
In this approach, the field operator # is represented as a 
sum of the (c- number) macroscopic component \&o, con- 
taining excitations with momenta k < ko -C ■f" 1 (where 
£ = h/y/mgn is the healing length) and a small opera- 
tor component 8^f describing excitations with larger mo- 
menta, i& = ^0 +8^. The momentum fco is chosen such 
that most of the particles are contained in H>q, however, 
its details do not enter into the lowest order corrections 
to g( 2 \r), which are 0(8^>) 2 . Using Wick's theorem, and 
the property of the thermal density matrix that (8^f) = 0, 
the pair correlation function is then reduced to 



,(2) 



(r) 



Re(<5* t (r)5*(0)) + Re(<5*(r)<5*(0)) 



(46) 



The normal and anomalous averages (Sijj^ (r)Sijj(O)} and 
(5i/j(r)5ip(0)) are calculated using the Bogoliubov trans- 
formation 



u k a k e 



v k a{e- ik * 



(47) 



where L is the length of the quantization box, a k and a\ 
are the annihilation and creation operators of elementary 
excitations, and (wfc,i>fc) are the expansion coefficients 
given by 



Uk 



efc — E k 



(48) 



and satisfying u\ 



1. Here e k = v / E k (E k + 2gn) is 



the Bogoliubov excitation energy, E k = h 2 k 2 /(2m), and 
we note that the following useful relationships between 
E k and e k hold: 



Ek 



\Je% + (gnf - gn, 



k 2 



k 2 + (2/C) 2 



1/2 



(49) 
(50) 



where £ = hj ^Jmgn is the healing length. The equilib- 
rium occupation numbers of the Bogoliubov excitations 



are given by n k 



(a{a k ) 



Applying the Bogoliubov transformation to the normal 
and anomalous averages in Eq. lj46|) gives 



g {2) (r) = 1 + — / dk cos(fcr) 



1 



x [(u k - v k ) 2 n k + v k (v k - Uk)] ■ (51) 

Using next Eq. lj48j) for the coefficients u k and v k we 
obtain the following result for the pair correlation func- 
tion 



1 



9 [Z} {r) = l + — 
Zirn 



dk 



— (2n fc + l)-l 



cos(fcr). 



(52) 

For convenience, we split the g( 2 )(r)-function into two 
parts corresponding to the contributions of thermal and 
vacuum fluctuations, 



<? (2) (r) = l + G (r) + G T (r), 



with 



Go(r) = — 
Zirn 



dk 



Ek 



— - 1 



cos(fcr), 



(53) 



(54) 



and 



+00 

Gr(r) = — [ dk—n k cos(kr). 
irn J e k 

—00 



(55) 
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We first evaluate the vacuum contribution Go(r), Eq. 
((54"]) . As shown in Appendix El the integral in (f54|) can 
be obtained exactly in terms of special functions, giving 



Go(r) = -V7[L-i(2V7w) - / 1 (2 v ^nr)] 



(56) 



where L_i(x) is the modified Struve function and h(x) 
is a Bessel function. The correlation length scale here is 
set by the healing length £ = h/^mgn = l/y^n. 



A. Quasi-condensate at low temperatures 

At very low temperatures when the excitations are 
dominated by vacuum fluctuations, whereas the thermal 
fluctuations are a small correction, the Gr(r)-term is cal- 
culated as follows. First, we substitute the explicit ex- 
pression for hk into Eq. l(55)l . giving 



+00 

G T (r) = — I dk^- 

irn I e h e e >" 1 — 1 



1 



cos(kr). (57) 



As shown in Appendix El for T -C gn (or t <C 7) the 
integral can be simplified and gives 



G T {r) 



2^7 



n 2 7r 2 r 2 47 



t 21 TTTnr\ 

— — cosech 



(58) 



Combining Eqs. (|53)l . l(56j) and ([58]) we obtain the 
following final result for this regime (t<7< 1): 



9 i2) {r) 



l-V7[L-i(2r/0--fi(2r/0] 

V7$ 2 _ . 2 / 7TTT 

27rr 2 8 7 3 /2 &ml v 2 7 £ 



(59) 



In the limit of r — > 0, the terms in the second line 
of Eq. (J39j) cancel each other and the large distance 
(r>?) asymptotics of the difference of special functions 
L_i(a;) — h( x ) ~ l/87ra; 2 ensures the expected inverse 
square decay of correlations^. At small but finite tem- 
peratures, the same large-distance asymptotics exactly 
cancels the inverse square behavior in the second line of 
Eq. lf59"j) leaving only the exponential decay 



9 {2 \r) 



8 7 3 / 2 



(60) 



to the uncorrelated value of g^ 2 \r) = 1. This is again in 
full agreement with the Luttinger liquid theory gjj. We 
note that even at T = 0, oscillating terms are absent, in 
contrast to the strongly interacting regime of Sec. IVI C[ 
Eq. [|71]). The limit r — > in Eq. ((59]) reproduces the 
result of Eq. (9) of Ref. Q, <?( 2) (0) = 1 - + 
ttt 2 /(24 7 3 / 2 ). In Fig.UJa) we plot Eq. {59]) for different 
values of the interaction parameter 7, and we note that 
the finite temperature correction term is negligible here. 




rK=1/ny 1/2 ] 
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FIG. 6: Nonlocal pair correlation g (r) in the weakly in- 
teracting regime, with r in units of the healing length £ = 
1/^/771: (a) low-temperature weekly interacting gas at r < 
7 <C 1, Eq. l[59|l : (b) weakly interacting gas at 7 <C r <C y 7 ^, 
Eq. 



B. Thermally excited quasi-condensate 

In the opposite limit, dominated by thermal rather 
than vacuum fluctuations and corresponding to 7 -C r <C 
y^, the thermal part of the pair correlation function is 
calculated as follows. We first note that large thermal 
fluctuations correspond to rife ^> 1, which in turn re- 
quires £fe/T <C 1. Thus, we replace rife in the integral 
(Ml by rife = [exp(efc/T) - l]" 1 ~ T/e k > 1. With 
this substitution, the integral for Gt{t) is dominated by 
the free-particle (quadratic in k) part of the Bogoliubov 
spectrum and the calculations in Appendix lO yield 



G T (r) 



2^7 



(61) 



This result is valid for r/£ < 1. For r/£ >■ 1 the main 
contribution to the integral in Eq. (j55]) comes from the 
phonon (linear in k) part of the Bogoliubov spectrum and 
one recovers the behavior given by Eq. (ffiTJJl . 

Combining Eqs. (|53]). (|56]) and 1(61]) we obtain the fol- 
lowing final result for this regime (7 <C r -C -^7 and 

r<0 = 



s (2)( r ) = i + .JL e -ar/« 



2^7" 

V7[L_x(2r/0-/i(2r/£)] 



(62) 
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The last two terms are due to vacuum fluctuations and 
are a negligible correction here, so the leading term gives 
an exponential decay of correlations [see Fig. EJb)] with 
a characteristic correlation length given by the healing 
length £ = Xj^Hn. The peak value at r = is g^(0) = 
1 + 7-/(2^/7), in agreement with Ref. [30]. 



VI. STRONGLY INTERACTING REGIME 

[7 > max{l, y 7 ?}] 

A. Perturbation theory in I/7 

By mapping the system onto that of a weakly attrac- 
tive ID fermion gas [6(J one can perform perturbation 
theory in I/7 <C 1. The formalism is the same as in Sec. 
IIV Al except that \& is now a fermionic field and the inter- 
action term in the Hamiltonian (J2j has to be modified to 
describe effective attractive interaction between fermions 
with matrix elements (in fc-space) V k = —2h 2 k 2 /(mnj) 
[63|. Then 

3 (2) M=9?2ooW + A<7 (2) M 

with 57=00 (r) = 1 — e~" Tr I 2 . The first order corrections 
to g^ir) are given by the Hartree-Fock approximation 
as a sum of the direct and exchange contributions 

A<?f (r) = fda J g V k T(k,a,r = 0)T(-k,a,r = 0)e ik \ 

(63) 



( 2 >M = - 



Ag^>(r) 



da J — V k T(k,a,r)T(-k,a,-r)e' 



kr 



(64) 



where 

T(k, a,r) = J dp G p+k (a)G p {-a)e lpr /27r, (65) 
in terms of the Green's function Gfc(cr) for free fermions. 

B. Regime of high-temperature "fermionization" 

We proceed with evaluation in the regime of high- 
temperature "fermionization" at temperatures well above 
quantum degeneracy, r 3> 1. In this regime, we use 
the Maxwell-Boltzmann distribution of quasi-momenta 
as the unperturbed state. In the temperature interval 
1<tC7 2 , the characteristic distance related to the in- 
teraction between the particles - the ID scattering length 
a\D = h 2 /mg ~ l\ja ~ 1/771 - is much smaller than the 
thermal de Broglie wavelength At, and the small pertur- 
bation parameter is am /At -C 1 [301 ] . 

>From the same formalism as in Sec. IIV A[ the free 
fermion Green's function is now given by 



G k (a) 



exp[(/3 + cr)(/i-ft 2 fc 2 /2m)], 
— exp[/ier — crfi 2 fc 2 /2m], 



-f3 < a < 0, 
< a < f3, 
(66) 



so the integral for T(k, a, r), Eq. {65]), gives 

r(fe a r) = _ ne - CT (^-^)' i2fe2 /2m / 3 e -mr 2 /(2/i 2 /3) e -ifcr ( T//3 

(67) 

Substituting Eq. (JBTJ) into Eqs. <[63j) and {64]) we obtain 
(see Appendix [Pi) 



A<?f (r) 



2Tn|r| 

7 



2 2 In 

-n tt J A 



-Mr), 



W) 



-Mr), 



(68) 
(69) 



(2) 

The only effect of the exchange contribution Ag\ is to 
cancel the delta-function in the direct contribution. This 
leaves us with the following result for the pair correlation 
function in the regime of high-temperature fermioniza- 
tion (1«t« 7 2 ): 



gV>(r) = l 



TTT 

U/ 7M 



^(rV27r/A T ) 



(70) 



In the limit r — > this leads to perfect antibunch- 
ing, ff( 2 )(0) = 0, while the small finite corrections (as 
in Ref. [30], g^ 2 \0) = 2r/ r ) 2 ) are reproduced at order 
7 -2 . The correlation length associated with the Gaus- 
sian decay of correlations in Eq. ((70]) is given by ther- 
mal de Broglie wavelength At = y/ in / (rn 2 ) . For not 
very large 7, the correlations do not decay in a sim- 
ple way, but instead show an anomalous, non-monotonic 
behavior with a global maximum at at r max — r )/2m. 
This originates from the effective Pauli-like blocking at 
short range and thermal bunching [g^ 2 '(r) > 1] at long 
range. As 7 is increased the position of the maximum 
diverges and its value approaches 1 in a non-analytical 
way gW(r max ) ~ 1 + (4t/ 7 2 ) exp(-7 2 /8r). 



Figure [2a) shows a plot of Eq. {70]) for various ra- 
tios of 7 2 /t. For a well-pronounced global maximum, 
moderate values of 7 2 /r are required (such as 7 2 /t ~ 5, 
with t = 8, 7 = 6), and these lie near the boundary 
of validity (7 2 /t 3> 1) for our perturbative result in the 
high-temperature fermionization regime. Exact numeri- 
cal calculations described in Ref. [34[ , and in more detail 
below in Sec. I VIII do. however, show qualitatively similar 
global maxima. 



C. Zero- and low-temperature (Tonks-Girardeau) 
regime 

At T = the procedure is straightforward and 
yields the known 0, H3] result 

(2) sin 2 (Q 4 sin 2 (Q 2tt 8 sin 2 (Q 



2 8_ 

' 7 d( 



C 

sin(C) 

c 



7 C 2 



dt sin(£i) In 
1 1—t 



7 dc e 

l + tl 



(71) 
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FIG. 7: Nonlocal pair correlation g (2 \r) as a function of the 
relative distance r in the strongly interacting regime, 7> 1: 
(a) regime of high-temperature "fermionization", 1 « t « 7 2 , 
Eq. H70p . with r in units of the thermal de Broglie wavelength 
At = \/47r/ (tti 2 ); (b) low temperature Tonks-Girardeau 
regime, Eq. I|7ip . for r = 0.01, with r in units of mean in- 
terparticle separation 1/n 




l I I T 
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r[A T ] 

FIG. 8: Behavior on the verge of the high-T fermionization 
regime for "f 2 /r = 4. The dashed line is Eq. (I70I) . 



in the gas. This can be interpreted as a quasi-crystalline 
order (with a period of ~ 1/n) in the two-particle sector 
of the many-body wave function even though the density 
of the gas is uniform. 

The oscillatory behavior of the pair correlation in this 
strongly interacting regime is similar to Friedel oscilla- 
tions in the density profile of a ID interacting electron 
gas with an impurity [HI]. We also mention that our 
derivation of Eq. ([71]) is equally valid for strong attractive 
interactions, i.e., when 7 < and I7I ^> 1, and therefore 
it describes the pair correlations in a metastable state 
known as super- Tonks gas (62I ]. 



D. Numerical results 



where C = imr. The last term here diverges logarithmi- 
cally with C and can be regarded as a first order pertur- 
bation correction to the fermionic inverse square power 
law. Accordingly, Eq. (|7Tj) is valid for ( <C cxp(7). 

At temperatures well below quantum degeneracy, r <C 
1, finite temperature corrections to Eq. ([7TJ are obtained 
using a Sommerfeld expansion around the zero temper- 
ature Fermi-Dirac distribution for the quasi-momenta. 
For rn <C t -1 this gives an additional contribution of 
t 2 sin 2 (7mr)/127r 2 to the right hand side of Eq. ([7l]) . 
which is negligible compared to the T = result as r <C 1 . 
At r = 0, Eq. ([TTjl gives perfect antibunching g^(0) = 0, 
which corresponds to a fully "fermionized" ID Bose gas, 
where the strong inter-atomic repulsion mimics the Pauli 
exclusion principle for intrinsic fermions. By extending 
the perturbation theory to include terms of order 7 -2 we 
can reproduce the known result for the local pai r corre- 
lation at zero temperature g (2) (0) = 47r 2 /37 2 [29l.[3oT|. 

In Fig. [7l[b) we plot the function gW(r), Eq. (JTTJ) , for 
various 7. According to the physical interpretation of the 
pair correlation function g^ 2 '(r), its oscillatory structure, 
and hence the existence of local maxima and minima at 
certain finite values of r, implies that there exist more 
and less likely separations between the pairs of particles 



Numerical calculations with the gauge- P method are 
able to reach only the low-7 (or, equivalently, high t) 
edge of the high-temperature fermionization regime, how- 
ever a comparison with Eq. 170]) is instructive. In Fig. [8] 
we see that the length scale on which antibunching occurs 
is still qualitatively given by Eq. I|70p while any discrep- 
ancies are of the same size as at r = 0. This is actually a 
general feature in all the parameter regimes explored by 
the numerics. Overall, the discrepancy between the I/7 
perturbation expansions lj39j) . |45j) . |70j) . and the exact 
behavior of g^ 2 '(r) at nonzero r is roughly the same as at 
r = 0. Since a calculation of <7^(0) [30] from the exact 
solution of the Yang- Yang integral equations is usually 
more straightforward to evaluate than the full stochastic 
calculation of g^ 2 \r), it can serve as a useful guide to 
whether a numerical calculation is warranted or not. 



VII. CLASSIC AL/FERMIONIZATION 
TRANSITION AND CORRELATION MAXIMA 

Figure [9] shows the behavior in the transition region 
between the decoherent classical and high temperature 
fermionization regimes (found with the gauge-P numeri- 
cal method), when one is far above the degeneracy tem- 
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FIG. 9: Crossover from decoherent classical to high tempera- 
ture fermionization regimes at high temperature. 
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FIG. 10: Situation when T > Td and the local second-order 
coherence is apparently unity. All curves plotted correspond 
to parameter values for which g' 2 '(0) = 1 in the crossover 
region between the classical decoherent and high-temperature 
fermionized gas. The dots (rather than triple lines here, for 
clarity) indicate la error bars. 

perature Td. One sees the appearance of a maximum 
in the correlations at finite range as the transition is ap- 
proached. As pointed out in Sec. I VI B| this arises from an 
interplay of thermal bunching and repulsive antibunching 
on comparable scales. A comparison of relevant length 
scales indicates that the t w here corresponds to 
At ~ aiD, where am is the "ID scattering length" that 
describes the asymptotic behavior of the wave function 
in two-body scattering. 

An interesting behavior occurs in the crossover regime 
when 7 2 /r ~ 0.1 — 0.4. Here we can have <r 2 '(0) = 1 just 
like in the quasi-condensate or "Gross-Pitaevskii" regime, 
indicating local second-order coherence. However, unlike 
the quasi-condensate regime, the non-local correlations 
on length scales of ~ A^ are not coherent, and in fact ap- 
preciably bunched. This is shown in Fig.[l0l It is a symp- 
tom of the broader correlation maximum phenomenon. 

The height of this maximum for more general parame- 
ters is shown in Fig. [TT]as a function of both g^{0) and 
7 2 /t. One sees that this behavior is well pronounced in 
the crossover between high temperature fermionization 
and decoherent classical regimes, peaking when g^ 2 '{0) — 
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FIG. 11: Heights of the anomalous peak of g^ 2 '(r) that 
occurs at nonzero r max , for different values of r, as func- 
tions of g^ 2 \0) - (a) and 7 2 /t - (b). The height is taken 
to be h = j' 2 '(r m „) — g^ 2 '(Q) at high temperatures when 
g (2) (0) > 1, and h = ff (2) (r max ) - 1 when ff (2) (0) < 1. The 
two regimes are separated by the dot-dashed vertical line in 
(a). Analytic results from Eq. lf39)l in the decoherent quantum 
regime are shown as a dashed line. Dots (rather than triple 
lines here, for clarity) indicate la error bars on the numerical 
results. 

1 (a situation shown also in Fig. [10]), or, equivalently, 
7 2 ~ 0.3r. As one reaches degenerate temperatures, 
the maximum peak height is reduced, and presumably 
disappears completely by the time the quasi-condensate 
regime is reached by going to smaller values of 7. Al- 
though we were unable to numerically reach the relevant 
quasi-condensate region for r < 1, a more refined numer- 
ical setup that improves the importance sampling or the 
/i(T) trajectory described in Appendix [A] may allow this. 



VIII. NUMERICAL LIMITATIONS 

Figure [32] shows the regime that was accessible using 
the relatively straightforward numerical scheme that was 
employed here, and detailed in Appendix [A] (It is the 
region above and to the left of the asterisks). In par- 
ticular, one sees that of the physical regimes described 
in previous sections, the decoherent classical, as well as 
parts of the decoherent quantum and high-temperature 
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FIG. 12: Regimes and their numerical accessibility: the as- 
terisks indicate the lowest r and highest 7 reachable using 
the gauge- P method as described in Appendix [A] The dark 
dashed line indicates the point at which (/ 2 '(0) = 1. 



fermionization regimes were accessible, while the quasi- 
condensate and Tonks-Girardeau regimes were not. 

The principal difficulty that is encountered, generally 
speaking, is the growth of statistical noise with increasing 
/?, i.e. decreasing r, which eventually prevents one from 
obtaining values of g^ (r) with a useful resolution. This 
arises in two different ways depending on the region of 
interest. 

Firstly, in the strongly interacting (fermionized) re- 
gion, one needs a correspondingly large coupling constant 
g oc 7 which leads to a relative increase of the impor- 
tance of the noise terms of the da.j /df3 equations in 

(|2T|l . This leads to large statistical uncertainty in the 
themselves or to the weight fl whose evolution depends 
on them. The upshot is that the inverse temperature 
(3 at which the noise becomes unmanageable becomes 
smaller and smaller as 7 grows. Technical improvements 
are unlikely to make a large dent in the problem in the 
fermionized regime because it ultimately stems from the 
fact that coherent states are no longer a good basis over 
which to expand the density matrix. They are not close 
to the preferred eigenstates of the system. Instead, one 
can think of constructing a phase-space distribution that 
uses a non-coherent-state basis, for example, a Gaussian 
basis [63| . This general approach - together with symme- 
try projections - has been utilized in successfully calcu- 
lating ground state properties of the strongly correlated 
fermionic Hubbard model (64j . 

Secondly, in the low 7 and r region, one has a different 
underlying source of statistical uncertainty. The longest 
relevant length here is either the coherence length 1$ or 
the healing length £, and for correct calculations in the 
large uniform gas one must simulate a system of a total 
size appreciably greater than these lengths. This in turn 



imposes a minimal total particle number 

N > max [C(4/t), 0{2/^j)}. (72) 

The thermal initial conditions of Eq. lfl~2j) lead to varia- 
tion in TV among trajectories, and since the Gibbs factor 
K (see Eq. lfT6|) ) grows linearly or faster with N, one 
also obtains a growing variation of K (y). This enters the 
dfl of Eq. lj2T]l and leads to a spread of the weights Cl(t) 
that grows rapidly (note the exponential growth of ft) 
with increasing N. However because of the long length 
scales, via ((72)1 . large N is needed to make accurate calcu- 
lations when t or 7 are much smaller than one. The end 
result is domination of the whole calculation by one or a 
few trajectories with the highest weight, for all realistic 
ensemble sizes S. 

As a corollary, significantly lower temperatures, even 
down to the quasi-condensate regime, are accessible at 
small 7 if one is prepared to sacrifice the assumption of 
an infinite-sized gas and consider periodic boundary con- 
ditions on some length L that is smaller than or compara- 
ble to the coherence/healing lengths. This approach was 
taken, e.g., in j6^|. This stops the rise of overall particle 
number, hence one has a much smaller spread of Gibbs 
factors Q among the trajectories, and in the final analysis 
- reduced statistical uncertainty. Such calculations are 
no longer as general, though, and are not considered in 
this paper. 

We would like to point out that the limitation in 
this regime may be overcome or alleviated if the rather 
simplistic importance sampling used in the numerical 
method were to be improved. The leading candidate is an 
improved importance sampling algorithm, possibly using 
a Metropolis sampling procedure, as outlined at the end 
of Appendix IA3I 

Finally, it is also possible that a more refined choice of 
(considered in Appendix IA 5ft may lead to somewhat 
improved coverage of the parameter space in general. 

IX. OVERVIEW AND CONCLUSION 

In conclusion, we have surveyed the behavior of the 
spatial two-particle correlation function in a repulsive 
uniform ID Bose gas. We have analyzed numerically the 
pair correlation functions for all relevant length scales, 
with the exception of several low-temperature transition 
regions (see Fig. [12] below the asterisks) which were not 
accessible by the numerical scheme we employed. Ap- 
proximate analytic results and methods have been pre- 
sented for parameters deep within all the major physical 
regimes. The key features of this behavior include: 

• Thermal bunching with <7^(0) ~ 2 and Gaussian 
drop-off at ranges At in the classical decoherent 
regime. 

• Exponential drop-off of correlations from </ 2 **(0) — 
2 at ranges 1$ in the decoherent quantum regime, 
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along with Gaussian-like rounding at shorter ranges 
~ At- 

• Suppressed density fluctuations with g^ 2 '{0) ~ 1 
and exponential decay at ranges of the healing 
length £ in the quasi-condensate regime. 

• Antibunching with g^(0) < 1 and Gaussian decay 
at ranges A T in the high-temperature fermioniza- 
tion regime. 

• Antibunching with g^(0) < 1 and oscillatory de- 
cay on ranges of the mean interparticle separation 
1/n in the Tonks-Girardeau regime. 

• Bunching at a range of ~ 0.3 At in the crossover 
between classical and fermionized regimes around 



7 



0.3t. 



Let us consider the regimes in turn, starting from 
the classical decoherent gas, then going anti-clockwise 
in Fig. [I2l The classical decoherent gas is well approxi- 
mated by Boltzmann statistics and is dominated by ther- 
mal fluctuations. The pair correlation function shows 
typical thermal bunching and a Gaussian decay, with the 
correlation length given by the thermal de Broglie wave- 
length At- 

As one reduces the temperature, the gas becomes de- 
generate, the thermal de Broglie wavelength becomes 
larger than the mean interparticle separation and loses 
its relevance. The correlation length increases and one 
enters into the decoherent quantum regime. Here, the 
dominant behavior of the gas is the ideal Bose gas bunch- 
ing, g( 2 \0) ~ 2, with large density fluctuations that de- 
cay exponentially on the length scale given by the phase 
coherence length l<p. Notably, the exponential behaviour 
starts to appear well above degeneracy first in the long- 
distance tails, being visible even around t ~ 10 as in 
Fig.H 

Reducing the temperature even further, while still at 
7 « 1, one enters into the quasi-condensate regime, in 
which the density fluctuations become suppressed and 
g( 2 )(0) ~ 1. In the hotter sub-regime dominated by ther- 
mal fluctuations, the pair correlation shows weak bunch- 
ing, <?( 2 '(0) > 1, while in the colder sub- regime domi- 
nated by quantum fluctuations one has weak antibunch- 
ing, g( 2 \0) < 1. In both cases the pair correlation decays 
on the length scale of the healing length £. 

We now move to the right on Fig. [TJl into the regime 
of strong interactions, while staying at temperatures well 
below quantum degeneracy, r <C 1. This is the Tonks- 
Girardeau regime, in which the density fluctuations get 
further suppressed due to strong interparticle repulsion. 
Antibunching increases and one approaches </ 2 '(0) = 
due to fermionization. The only relevant length scale 
here is the mean interparticle separation, 1/n, and the 
pair correlation function decays on this length scale with 
some oscillations. 

We next move up on Fig. [TJl to higher temperatures, 
and enter the regime of high-temperature fermionization. 



At short range, the pair correlation here is still anti- 
bunched due to strong interparticle repulsion, however, 
thermal effects start to show up on the length scale of At- 
As a result of these competing effects, the nonlocal pair 
correlation develops an anomalous peak, corresponding 
to bunching at-a-distance, with g^ 2 ^(r max ) > 1, begin- 
ning around r ~ 7 2 /2. 

As we increase the temperature even further, the ther- 
mal effects start to dominate over interactions and the 
antibunching dip gradually disappears. At temperatures 
r ~ 7 2 we observe a crossover back to the classical deco- 
herent regime. 

Our results provide new insights into the fundamental 
understanding of the ID Bose gas model through many- 
body correlations. Calculation of these non-local correla- 
tions is not accessible yet through the exact Bethe ansatz 
solutions. We expect that our theoretical predictions will 
serve as guidelines for future experiments aimed at the 
measurement of nonlocal pair correlations in quasi-lD 
Bose gases. 
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APPENDIX A: TECHNICAL APPENDIX FOR 
THE GAUGE-P CALCULATIONS 

1. Instability of the stochastic equations and its 
removal with a stochastic gauge 

A straightforward application of the ungauged diffu- 
sion Eqs. lfT9|) is foiled by the presence of an instabil- 
ity in the daij /df3 equations. We can see this if we 
first consider the evolution of Nj and discard the noise 
and kinetic-energy parts of the equation. Taking the de- 
terministic part from the Stratonovich calculus which is 
used for our numerics (this introduces the 1/2 term be- 
low), one has 



dp 



N 4 



Ax 



Nj — — 
J 2 



(Al) 



There are stationary points at the vacuum Nj = and 
at Nj = N a = 1/2 + [i e lS.xlg, with the more positive sta- 
tionary point (usually N a ) being an attractor, and the 
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more negative a repellor [see Fig. [13] (a) ]. The determin- 
istic evolution is easily solved, and starting from a time 
(3a gives later evolution as 



(A2) 



If has a negative Nj(Po), which is possible due to the 
action of the noises C, then at a later time 
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sing 
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Me 
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(A3) 



the solution has diverged to negative infinity. This be- 
havior of the deterministic part of the equations is known 
as a "moving singularity" and is a well-known indicator 
of non- vanishing boundary terms when an integration-by- 
parts is performed on the operator equation (fl4|) [isLIBa]. 
It implies that the FPE (fl8l) is not fully equivalent to 
quantum mechanics. 

The use of a stochastic gauge to remove this kind of 
instability has been described in [47], and in more detail 
in [!||. The gauge identity, Eq. ifTTj) . can be used on Eq. 
(fT4"ll to introduce an arbitrary modification to the de- 
terministic evolution (arising from first order derivative 
terms) for the price of additional diffusion in the weight 
ft. Since the gauge identity is zero, we can add an ar- 
bitrary multiple of it to Eq. (fl"4"l) . In particular, if we 
add 




with arbitrary functions Qj{v,f3), and perform the sub- 
sequent steps as before, then the diffusion matrix in the 
resulting FPE remains positive semidefinite (no negative 
eigenvalues), and the resulting Ito diffusion equations of 
the samples become 
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instead of fl9|). The atj equations are modified and com- 
pensating correlated noises have been added to the Q 
equation. 

We now wish to choose the functions Gj , called stochas- 
tic gauges, so that the instability is removed, keeping 
also in mind the goal of keeping the (now unbiased) sta- 
tistical uncertainty manageable. Heuristic guidelines for 
choosing gauges have been investigated in detail in [4§1 |. 
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FIG. 13: Deterministic phase space for Stratonovich for of the 
dNj equation, when /i e = 0. (a): ungauged, (b): using the 
gauge (j A6[) . The moving singularity in (a) is shown with a 
large arrow, the attractor in (b) at \Nj\ = N a with a thick 
dashed line. 



Several choices for a single-mode system were also inves- 
tigated there in Sec. 9.2 in terms of resulting statistical 
uncertainties. The aim is to remove the real part of Nj 
from the atj equation when it is negative, so as to neutral- 
ize the moving singularity. While for a single mode the 
"radial" gauge was found to give the best performance, 
later tests that we have performed on the full multimode 
(M ^> 1) ID gas show that the "minimal" drift gauge 



^ = i(ReJV,--|^.|) 



9 



2Ax 



(A6) 



gives better performance for this system. This is because 
it introduces the smallest modifications needed to remove 
the moving singularity, and hence the smallest noise con- 
tributions to the weight fl. The weight becomes much 
more important for multimode systems because each of 
the M modes adds its own contribution to it, the total 
of which can become large. The phase-space modifica- 
tion for a single mode for the ungauged Eq. I|A1|) and 
gauged equations is shown in Fig. [131 One sees that in 
the "classical" Re[ATj] 3> Im[iVj] region the trajectories 
are practically unchanged. The final Ito equations to be 
integrated are i[2Tj) . Comparisons to known exact results 
such as energy and density [3], and </ 2 ^(0) indicate 
no deviations beyond what is predicted by the unbiased 
statistical uncertainties, Eq. l(26|) . with the new gauged 
equations. Such a comparison can be seen in Fig. 2 of 
Ref. H. 



2. Integration procedure 

The actual integration is performed using a split-step 
semi-implicit method described in [66j|, which requires 
the use of the Stratonovich stochastic calculus. There, it 
was shown to be highly superior to other low-order meth- 
ods in terms of stability. Although a low order Newton- 
like method, with the right choice of variables its per- 
formance is remarkably good. High-order methods such 
as Runge-Kutta or others suffer from serious complica- 
tions when noise is present. In particular, one has to be 
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very meticulous in tracking down and compensating for 
all the non-zero correlations within a single time-step — 
these are much more complicated than the simplest cor- 
rection terms appearing in the Stratonovich semi-implicit 
method used here. 

Due to the multiplicative form of the equations (|2T|) . 
it is highly advantageous to use logarithmic variables, 
which is made possible if one uses a split-step method. 
Here, a A/3 timestep consists of the following four stages: 
First the interaction part (containing g) is integrated in 
real space over a time-step A/3. Second, the fields are 
Fourier-transformed to fc-space, giving a^ v '(k). Thirdly 
the kinetic-energy contributions are integrated over A/?, 
and finally one Fourier-transforms back into real space, 
ready to start the next timestep. The Stratonovich 
gauged evolution equations for the real space stage are 



a universal common factor in the / dv integral. Hence, if 
we manually scale the weight ft by some factor F(v) of 
our choice: ft — > Q'F(v), and simultaneously rescale the 
distribution according to G(v) — ► G \v)j F(v), then with 
ft' and G' one obtains exactly the same results in the 
infinite-number-of-samples limit as with Gft. However, 
the actual samples are differently distributed, which is 
advantageous for finite sample numbers. To reduce the 
weight sampling problem, one wants to make such a mod- 
ification F(v) that both G'(v) and fl'G'(v) peak in the 
same region of the phase space of v. 

To proceed, it is convenient to consider Fourier- 
transformed variables in fc-space, where the non- 
interacting evolution can be easily exactly solved. Define 
then 
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where k takes on discrete values from —ir/Ax to tt/Ax. 
The "naive" initial distribution (TT2l) then becomes 
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E(Re^-|iV,|)cf(/3) (A7a) G (v) = 6 2 (lnQ) ]j6 2 (a k («+)* )- 
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while for the fc-space stage they are 



This is a thermal distribution which is uniform over all 
k. The ideal gas (i.e. g = 0) evolution of equations l|2"Tj) 
then leads to 
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where 



3. Importance sampling 



5 fc (0) 



n x Tjk, 



(All) 



The simulated equations lj2lj) include evolution of both 
the amplitudes cq and weight Q. This combination can 
cause sampling problems for observable estimations, Eq. 
(|22|l . when maximum weights occur for very rare trajec- 
tories. As it turns out, this was a serious issue for the 
majority of calculations reported here because while the 
initial distribution (fl2|) samples the P = system well, 
this is not necessarily the case during the later evolution 
into P 3> that is of most interest. Fortunately, fairly 
rudimentary importance sampling was able to deal with 
this for a wide range of parameters. 

The essence of this approach is to pre-weight trajec- 
tories in such a way that the part of the distribution 
with maximum weight fl coincides with the majority of 
samples at the target time of interest Pt, rather than at 
P = 0. The price paid is that the j3 = distribution is 
then poorly sampled, but this is not important to us as 
we are interested rather in the target P t - 

Pre-weighting is made possible because in all observ- 
able calculations l(22|). the combination [G(w)ft] occurs as 



with rjk being independent complex Gaussian noises with 
variance unity, (rf^rjk^s = $kk' ■ One can see that (|A10j) 
is not necessarily anywhere near a well-sampled ideal 
gas Bose-Einstein distribution at temperature p, which 
would have 
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with 



n f(p) = {exp [-[i(P)P + h 2 k 2 p/2m] - 1} 1 

being the usual Bose-Einstein distribution. 

For the purpose of the simulations presented here, 
a fairly crude yet effective importance sampling was 
applied as follows. For relatively weak coupling g, a 
very rough but useful estimate of the thermal state at 
coarse resolution is that the Fourier modes are decoupled 
and thermally distributed with some mean occupations 
n-k(Pt) at the target time Pt that we are interested in. In 
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practice we will choose some estimate of the guiding den- 
sity n k (p t ). The desired equal weight sampling at time 
Pt would then correspond to the distribution 



[CXk 



exp[-\a k \ 2 /n k (f3t)] 



(A13) 



which leads to samples given by a k = ^/rifc(/3) r\ k and 
= 1. What we are interested in is the correspond- 
ing distribution of samples at j3 = 0. An estimate of 
the initial distribution that leads to G est (v,Pt) can be 
obtained by evolving l|A13|) back in imaginary time us- 
ing only kinetic interactions. This is again rather rough, 
since deterministic interaction terms oc g are omitted, 
not to mention noise, but it is simple to carry out and 
proved sufficient for our purposes here. One obtains then 
an estimated sampling distribution for samples at (3 = 0: 



exp(-|5,| 2 K amP ) 
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and the pre-weight Oo = 0(0) now depends on the set 
of particular values of a k at p = obtained for a given 
sample, according to 
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For most of the simulations reported here, taking 
rik(Pt) to be just the ideal gas Bose-Einstein distribu- 
tion n k d (Pt) was sufficient. However, once the chemical 
potential /i(/?t) approaches or exceeds zero, this estimate 
is no longer useful. A better choice for n k (f3t) is the den- 
sity of states function p k of the exact Yang and Yang 
solution 0, although it should be noted that this is not 
the density of actual particles that we seek. In practice, 
our approach was to first run a calculation based on this 
estimate n k {p t ) = Pk{Pt), obtain a better estimate of the 
real density from this full stochastic calculation by evalu- 
ating the expectation value of ^ k ^k using Eq. lj22j) . then 
finally use this expectation value to choose an improved 
preweighting function n k (/3 t ) for a "second-generation" 
calculation. 

One important point to make regarding the choice of 
nk(Pt) is that one should endeavor always to choose the 
preweighting guide density n k ((3 t ) equal or greater than 
the real density, never smaller. The reasoning behind 
this is as follows: Suppose first one chooses a n k (f3 t ) 
guiding function that is much smaller than the true k- 
space density n). me (/3t). This means that the variance of 



the a k samples will be too small to recover the physical 
value of the density upon averaging (\a k \ 2 tt)s without 
resorting to very large weights for the largest \a k \ sam- 
ples. In practice, if the ratio n k /n k IU ° is small, then the 
typical trade-off that occurs is that the largest contri- 
bution to O | 2 comes from those \a k \ that are many 
standard deviations from the mean. Their rarity is com- 
pensated for by a very large O. However, this is fatal for 
practical numbers of samples because in fact not even 
one of the samples one obtains ends up in this highest- 
contribution region at many standard deviations from 
the mean. For n k /n t k uc < 1/2, the number of samples 
with \a k \ 2 > n k will be oc S f[ k exp [-(n\ mc /n k ) 2 /2) , 
i.e. vanishing, leading to a systematic error. 

In contrast, the opposite situation when n k {p t ) is cho- 
sen too large is much more benign. Following the above 
reasoning, one gets a distribution of a k samples that is 
too broad, with the result that a majority of samples are 
too far away from physical values of |5fc| 2 and their exces- 
sive abundance must be compensated for by giving them 
a correspondingly small weight. However, for reasonably 
large numbers of trajectories, there always remains a core 
of the smallest samples that are in the region of most im- 
portant contributions. The number of these samples is 
of the order of SY\ k n k Ille /n k ((3 t ), which is reasonable in 
practice as long as the estimate n k (P t ) is not extremely 
poor. 

Finally, it should be mentioned that superior impor- 
tance sampling schemes to the crude one we have em- 
ployed here could be implemented and may allow one 
to reach much lower temperatures than presented here. 
A first step would be to keep the (3 = fit distribution 
estimate, Eq. 1|A13|) . but estimate the resulting initial 
samples at f3 = in a more accurate manner. To do 
this, one could choose the (3 = fit samples according to 
(Pt) = y/rik((3t) f]k and lnO(/3 f ) = as usual, but 
then evolve them back in time to j3 = numerically, 
using the deterministic part of the full equations (|2"T|) . 
This would give a superior estimate of the initial distri- 
bution as it takes into account g ^ mean field effects 
as well as kinetic evolution. Having these (3 = sam- 
ples, one would then proceed forward in time with the 
full stochastic evolution. 

A further refinement would be to choose initial (3 = 
samples via the Metropolis algorithm, so that the ini- 
tial samples v are distributed according to T \v\ , where 
T = |0(/3j)| when 0(/3 t ) is calculated according to the 
deterministic part of the evolution, Eq. (|2Tj) . starting 
from O(0) = 0. This avoids the arbitrariness of the 
crude Gaussian choice, Eq. (|A13[i . A final, but numer- 
ically intensive approach would be to sample the phase- 
space variables a.j{(3t) and Im[lnO](/3t) directly via a 
Monte Carlo Metropolis algorithm whose free parame- 
ters to be varied include both the initial noises rj k and all 
the time-dependent noises Q (P) for a given time lattice 

pe(0,p t ). 



19 



4. Trust indicators for sampling 

One should mention two heuristic trust indicators that 
we use extensively to exclude bad sampling of the under- 
lying phase-space distribution. 

Firstly, let us point out that the behavior of the evolu- 
tion equations 1|A7[I is such that one builds up an approxi- 
mately Gaussian distribution of the logarithmic variables 
(leaving aside the evolution of Nj itself, which is initially 
small). This means that the stochastic averages to be 
evaluated, e.g., in Eq. ([25]) . involve means of exponen- 
tials of approximately Gaussian random variables (as per 
to = (e u ) with v Gaussian). A feature of such means is 
that if the variance of the logarithm Re[w] exceeds a value 
of around 10 the mean m begins to have systematic error 
when calculated with a ny p ractical sample sizes. This 
is discussed in detail in |48l. f67| . As a result, when cal- 
culating observables with some expression (F(v))s, one 
must also check that the variance of its logarithm is small 
enough, i.e. that 



Vr = ((ln\F(v)\f) s -(\n\F(v)\)l < 10. 



(A17) 



If this is not satisfied, the results for (F(v))s must be 
considered suspect. 

Secondly, sampling problems of this sort usually make 
themselves visible if one compares two calculations with 
widely different sample sizes. In practice one can evalu- 
ate an average and its uncertainty with S samples, and 
with 5/10 samples (where, of course, 5/10 ^> 1). If the 
difference is statistically significant the result of the S 
sample average again should be considered suspect. 



5. Choice of intermediate n(f3) 

If one is primarily interested in the behavior of the 
system around some target temperature fa and chemical 
potential n(fa) (alternatively - density), then the values 
of ^i(fa) at intermediate times f3 < fit can in principle be 
chosen at will. 

In practice, however, some choices lead to smaller sta- 
tistical uncertainty than others because the intermedi- 
ate values of density affect the amount of noise gener- 
ated during the evolution. A preliminary investigation of 
fi choice in indicated some heuristic guidelines that 
were also followed in the present work: 

(i) It is advantageous to not vary /J. e (/3) too much over 
the course of the simulation. Excessive variation leads to 
increased noise. 

(ii) A constant or piecewise-constant value of /i e is also 
advantageous because the ideal-gas part of the evolution 
can then be calculated exactly in logarithmic variables 
(|A7bj) , and step-size is only important for the interaction 
part of the evolution. 

(Hi) It is advantageous to choose an initial density that 
is much smaller than the final one at fa both for sta- 
tistical sampling reasons and because this puts the ini- 
tial gas much further into the classical decoherent regime 



(t 3> 7 2 ), where the initial condition fT2|) applies, than 
the final regime. 

In practice, our simulations used the following form 



1 z(/3 + A/3) 



(A18) 



which is piecewise constant over a time step A/3, with 
the fugacity 



when (3 < fa, 
zt exp l-fff In fi , when (3 > fa. 

(A19) 

Here, fa and z t = e^ 1 " 3 ' are the target inverse temperature 
and fugacity, and fa and Zi are numerical constants for 
the initial high temperature state that we chose to be 
z 2 = z t 2 /1000 and fa = fa/WOO. 

Given the difficulty of precisely analyzing the statisti- 
cal behavior, it is unclear whether a wiser choice of /i(/3) 
may lead to significant improvements over the results pre- 
sented here. However, this is the most successful choice 
of those we tried. 



APPENDIX B: INTEGRALS IN PERTURBATION 
THEORY IN 7 

We begin with Eq. f32|) and substitute the expression 
for T(fc, er) in Eq. (37} to give 

r 2 m8 

ft V 7T J ^2/4-^-/3/2)2 



(Bl) 

Next we make the substitution t = (2/ fa](a — /3/2) and 
y = ryjmj '(ti 2 fa) to give 

A.9 (2) M = 



dt , (B2) 



9 rnj3 

n\ tt y_r" vi - ~£ 

9 I m P 



dx- 



l + x- 



; (B3) 



whe re the l ast equality follows from the substitution t = 
x/Vl + x 2 . The exponent in the integrand of Eq. (|B3j) 
can be represented as a Gaussian integral 



1 



dke -k +2ikyx_ 



(B4) 



Then, changing the order of integration in Eq. (|B3|1 we 
arrive at 



AgM(r) 




£ i2kyx 
1+X 2 



= - 2g ^ He-^dk. 



rdxdk 



(B5) 
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The final result shown in Eq. (|38|l follows trivially from 
a shift in the integration variable k — > k — \y\, and the 
definition of the complimentary error function, 



2 f°° 
erfc(|y|) = —= / dk 

V 1 J\y\ 



,-k 2 



(B6) 



APPENDIX C: INTEGRALS IN THE 
BOGOLIUBOV TREATMENT 

We first evaluate the vacuum contribution Go(r), Eq. 
(|54|) . Writing down the integral explicitly, in terms of k, 
and transforming to a new variable x = fc£/2, we have 



Go(r) 



2 





Integrating by parts, gives 

oo 

G (r) = — /dz 



- 1 



cos(2rx/0- (CI) 



3in(2^/7nra;) 

( 1+ x 2)3/2 ■ 



(C2) 



where we have introduced 7/ = gnx/T = e/T, 
Finally we make use of the following integral 

J2 



dy 



ycos(ay) 
ev - 1 



1 7T^ 

— cosech (ira) , (C8) 

2a 2 



and obtain Eq. (j58|) . 

In the opposite limit, dominated by thermal fluctua- 
tions and corresponding to 7 r 1, we first note that 
large thermal fluctuations correspond to fik ^> 1, which 
in turn requires e^/T <C 1. Thus, we replace fik in the 
integral JM]) by n fe = [exp(e fe /T) - l]" 1 ~ T/e fc > 1. As 
a result, the thermal contribution Gt{t) becomes 



G T (r) ~ — 
7rn 



+00 EfcT 

dk — — cos(fcr) 



4mT 
irh 2 n 



+00 



dk- 



cos(fcr) 

2 + (2/6 2 



h 2 n 



-2r/£ 



(C9) 



which is valid for r/£ < 1. Rewriting this in terms of the 
dimensionless parameters 7 and r we obtain Eq. (J6XJ) . 
For ?*/£ 3> 1 the cosine term becomes important and the 
values of momenta in the integral Eq. (|C4[) are cut off by 
1/r <C £. In this regime one can use the approximation 
that led to Eq. (fU8l) , 



The integral in <|C2[) can be expressed in terms of spe- 
cial functions [68], giving 



(C3) 



G (r) = -y/i [L^(2^nr) - h(2^nr)} 



The finite temperature term Gt{t), Eq. (J57j) , is 
evaluated by performing variable changes acc ording to 
E = h 2 k 2 /(2m), followed by e = y/E(E + gn) and then 
x = e/gn. In this way we transform the integral over k 
to an integral over x 



G T {r) 



2m g 
Tr 2 h 2 n 



dx 



VT 



1 



1 + x 2 



1/2 



cos[k(x)r] 

e gnx/T _ I ' 



(C4) 

where k(x) = [2mgn(Vl + x 2 — Yj/H 2 ] 1 / 2 . So far we have 
not made any additional assumptions or approximations. 

By inspecting the integrand in Eq. (|C4|) one can see 
that for T <C gn the main contribution to the integral 
comes from 1 <C 1. Therefore for T <C gn (r <C 7) we 
can simplify the integral by treating x in the integrand 
as a small parameter. Accordingly, we obtain 



VT 



1 



nl/2 



—=x, iCl, 
V2 



k{x) 



mgn 



h 2 



x, x <C 1, 



and therefore 



G T (r) 



47T7 3 / 2 



dy 



y cos(rnry /2^/j) 



ev - 1 



(C5) 
(C6) 

(C7) 



APPENDIX D: INTEGRALS IN PERTURBATION 
THEORY IN 1/7 

We begin by evaluating the direct contribution given 
by Eq. by substituting Eq. J67]) , 



da 



-1 IT 
7T7 V 2 



— 00 

1 



— J- / rfs 



dfc / 2ft 2 fc 2 

27T \ 771717 
'OO 



dqq e 



Zpiqy-sq (1-*) 



(Dl) 



where we hav e affected th e change of va riables g = /3s, 
g = ^//3h 2 /mk and 7/ = y/m/ \(3h 2 )r = yj (rn 2 /2)r. The 
integration with respect to q can then be done using in- 
tegration by parts, which yields 



(2) 



-1 It f 1 ,2s(l- s) ~y 2 ^_ y 2 /[4s{1 _ s)] 



ds 



4 7 V 2tt J q s 5 / 2 (1-s) 5 /2 



-1 /2r 



■fVv7.*l 1 -i%- 



V/(i-t 2 ) 



(1-i 2 ) 



3/2 



(D2) 

where the last equality follows from the substitution s = 
(t + l)/2. The simplest way to solve the integral in Eq. 
(p2|) is by comparison with Eq. (|B2[) in Appendix [Bj In 
doing so, one may observe 



dt 11- 



2y 2 



1 

2 



dy 



dt- 



1-t 2 
exp 



exp 



y 

i~t 2 



(l_ t 2)3/2 



(D3) 



v 



VT^f 2 



d 2 

dy 2 



7T— erfc(|y|). (D4) 
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The result shown in Eq. f68|) then follows trivially from 
this. 

In order to calculate the exchange contribution we be- 
gin with Eq. f64|) and substitute Eq. I|67p . which imme- 
diately yields 



where the first equality comes from the su bstitut ion s — 
(v + l)/2 and the second from v = t/yL + t 2 . Both 
terms are standard definite integrals it is straightforward 
to show that 



Agi 2 \r) = 



/2 F e (y / T n 2 r 2 /2) (D5) 



where F e (y) = J* ds J dqq 2 e~ s ^' s ^+^ 1 ~ 2s ^y /^f 2 , 
and s, q and y are defined the same was as for the direct 
contribution. The integration with respect to q can be 
carried out using integration by parts, leaving an integral 
with respect to s: 



A 5 f = -5{r). 



(D7) 



i exp 

ds 

3 sV 2 (l 

i exp 
dv 

-i 

oo 



4s(l-s) 



s )3/2 



2 2 

y_u 
1_„2 



J -V 2 ) 3 / 2 

dt [l - 2y 2 t 2 ] i 



2/ 2 (l-2 S ) 5 



2s(l 
2v 2 y 2 " 



(D6) 



Thus the only effect of the exchange contribution is to 
cancel the delta-function contribution coming from the 
direct contribution at r = 0. 
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